###Space, Time, Space vs. Time#####
###02/01/2022######################
###Generate Plots##################

#Generic working directory
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#Get packages, attach them
pks=c("ggplot2", "reshape2")
needed=setdiff(pks,installed.packages())

if( length(needed) > 0){
  options(Ncpus=parallel::detectCores() %/% 1.25)
  install.packages(needed)
}

sapply(pks, library, character.only=T)
session = sessionInfo()

###Spatially dominated X results
beta.results.static.all <- matrix(NA, nrow=7, ncol=6)
beta.results.ldv.all <- matrix(NA, nrow=7, ncol=6)
beta.results.sar.all <- matrix(NA, nrow=7, ncol=6)
beta.results.stadl.all <- matrix(NA, nrow=7, ncol=6)
phi.results.ldv.all <- matrix(NA, nrow=7, ncol=6)
phi.results.stadl.all <- matrix(NA, nrow=7, ncol=6)
rho.results.sar.all <- matrix(NA, nrow=7, ncol=6)
rho.results.stadl.all <- matrix(NA, nrow=7, ncol=6)
lm.test.all <- matrix(NA, nrow=7, ncol=6) 

load("SpatiallyDominatedXphi_0.RData")

beta.results.static.all[,1] <- beta.results.static[1,]
beta.results.ldv.all[,1]  <- beta.results.ldv[1,]
beta.results.sar.all[,1]  <- beta.results.sar[1,]
beta.results.stadl.all[,1]  <- beta.results.stadl[1,]
phi.results.ldv.all[,1]  <- phi.results.ldv[1,]
phi.results.stadl.all[,1]  <- phi.results.stadl[1,]
rho.results.sar.all[,1]  <- rho.results.sar[1,]
rho.results.stadl.all[,1]  <- rho.results.stadl[1,]
lm.test.all[,1] <- lm.test[1,]

load("SpatiallyDominatedXphi_0.1.RData")

beta.results.static.all[,2] <- beta.results.static[1,]
beta.results.ldv.all[,2]  <- beta.results.ldv[1,]
beta.results.sar.all[,2]  <- beta.results.sar[1,]
beta.results.stadl.all[,2]  <- beta.results.stadl[1,]
phi.results.ldv.all[,2]  <- phi.results.ldv[1,]
phi.results.stadl.all[,2]  <- phi.results.stadl[1,]
rho.results.sar.all[,2]  <- rho.results.sar[1,]
rho.results.stadl.all[,2]  <- rho.results.stadl[1,]
lm.test.all[,2] <- lm.test[1,]


load("SpatiallyDominatedXphi_0.2.RData")

beta.results.static.all[,3] <- beta.results.static[1,]
beta.results.ldv.all[,3]  <- beta.results.ldv[1,]
beta.results.sar.all[,3]  <- beta.results.sar[1,]
beta.results.stadl.all[,3]  <- beta.results.stadl[1,]
phi.results.ldv.all[,3]  <- phi.results.ldv[1,]
phi.results.stadl.all[,3]  <- phi.results.stadl[1,]
rho.results.sar.all[,3]  <- rho.results.sar[1,]
rho.results.stadl.all[,3]  <- rho.results.stadl[1,]
lm.test.all[,3] <- lm.test[1,]


load("SpatiallyDominatedXphi_0.3.RData")

beta.results.static.all[,4] <- beta.results.static[1,]
beta.results.ldv.all[,4]  <- beta.results.ldv[1,]
beta.results.sar.all[,4]  <- beta.results.sar[1,]
beta.results.stadl.all[,4]  <- beta.results.stadl[1,]
phi.results.ldv.all[,4]  <- phi.results.ldv[1,]
phi.results.stadl.all[,4]  <- phi.results.stadl[1,]
rho.results.sar.all[,4]  <- rho.results.sar[1,]
rho.results.stadl.all[,4]  <- rho.results.stadl[1,]
lm.test.all[,4] <- lm.test[1,]

load("SpatiallyDominatedXphi_0.4.RData")

beta.results.static.all[,5] <- beta.results.static[1,]
beta.results.ldv.all[,5]  <- beta.results.ldv[1,]
beta.results.sar.all[,5]  <- beta.results.sar[1,]
beta.results.stadl.all[,5]  <- beta.results.stadl[1,]
phi.results.ldv.all[,5]  <- phi.results.ldv[1,]
phi.results.stadl.all[,5]  <- phi.results.stadl[1,]
rho.results.sar.all[,5]  <- rho.results.sar[1,]
rho.results.stadl.all[,5]  <- rho.results.stadl[1,]
lm.test.all[,5] <- lm.test[1,]

load("SpatiallyDominatedXphi_0.5.RData")

beta.results.static.all[,6] <- beta.results.static[1,]
beta.results.ldv.all[,6]  <- beta.results.ldv[1,]
beta.results.sar.all[,6]  <- beta.results.sar[1,]
beta.results.stadl.all[,6]  <- beta.results.stadl[1,]
phi.results.ldv.all[,6]  <- phi.results.ldv[1,]
phi.results.stadl.all[,6]  <- phi.results.stadl[1,]
rho.results.sar.all[,6]  <- rho.results.sar[1,]
rho.results.stadl.all[,6]  <- rho.results.stadl[1,]
lm.test.all[,6] <- lm.test[1,]




##Figure A8a - LDV phi bias 

rho = seq(0.0, 0.30, 0.05)
phi.results.ldv.all = cbind(phi.results.ldv.all, rho)
colnames(phi.results.ldv.all ) <- c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50","rho")

cplot_dataA10 = reshape(as.data.frame(phi.results.ldv.all), 
                      varying = c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50"), 
                      timevar = "phi",
                      direction = "long",
                      idvar = "rho",  
                      sep="-"
)


cplot_dataA10$phi <- as.factor(cplot_dataA10$phi)

p = ggplot(data=cplot_dataA10, aes(x=rho, y=P, group=phi, linetype=phi)) 
p + geom_line(size = 0.75) +  
  theme_bw(base_size = 15) +
  theme(legend.key.width = unit(1.5,"cm")) +
  ylab(expression(hat(phi)~-phi))+
  xlab( expression(rho[Y])) +
  labs(linetype=expression(phi[Y]))

ggsave("FigureA8a.pdf", height=4, width=6, units="in", dpi=600)


##Figure A8b - LDV beta bias 
rho = seq(0.0, 0.30, 0.05)
beta.results.ldv.all = cbind(beta.results.ldv.all, rho)
colnames(beta.results.ldv.all ) <- c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50","rho")

cplot_dataA11 = reshape(as.data.frame(beta.results.ldv.all), 
                      varying = c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50"), 
                      timevar = "phi",
                      direction = "long",
                      idvar = "rho",  
                      sep="-"
)


cplot_dataA11$phi <- as.factor(cplot_dataA11$phi)

p = ggplot(data=cplot_dataA11, aes(x=rho, y=P, group=phi, linetype=phi)) 
p + geom_line(size = 0.75) +  
  theme_bw(base_size = 15) +
  theme(legend.key.width = unit(1.5,"cm")) +
  ylab(expression(hat(beta)[LDV]~-beta))+
  xlab( expression(rho[Y])) + 
  labs(linetype=expression(phi[Y]))

ggsave("FigureA8b.pdf", height=4, width=6, units="in", dpi=600)




##Figure A9a - SAR rho bias 


rho = seq(0.0, 0.30, 0.05)
rho.results.sar.all = cbind(rho.results.sar.all, rho)
colnames(rho.results.sar.all ) <- c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50","rho")

cplot_dataA12 = reshape(as.data.frame(rho.results.sar.all), 
                       varying = c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50"), 
                       timevar = "phi",
                       direction = "long",
                       idvar = "rho",  
                       sep="-"
)


cplot_dataA12$phi <- as.factor(cplot_dataA12$phi)

p = ggplot(data=cplot_dataA12, aes(x=rho, y=P, group=phi, linetype=phi)) 
p + geom_line(size = 0.75) +  
  theme_bw(base_size = 15) +
  theme(legend.key.width = unit(1.5,"cm")) +
  ylab(expression(hat(rho)~-rho))+
  xlab( expression(rho[Y])) + 
  labs(linetype=expression(phi[Y]))

ggsave("FigureA9a.pdf", height=4, width=6, units="in", dpi=600)


##Figure A9b - SAR beta bias 

rho = seq(0.0, 0.30, 0.05)
beta.results.sar.all = cbind(beta.results.sar.all, rho)
colnames(beta.results.sar.all ) <- c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50","rho")

cplot_dataA13 = reshape(as.data.frame(beta.results.sar.all), 
                       varying = c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50"), 
                       timevar = "phi",
                       direction = "long",
                       idvar = "rho",  
                       sep="-"
)


cplot_dataA13$phi <- as.factor(cplot_dataA13$phi)

p = ggplot(data=cplot_dataA13, aes(x=rho, y=P, group=phi, linetype=phi)) 
p + geom_line(size = 0.75) +  
  theme_bw(base_size = 15) +
  theme(legend.key.width = unit(1.5,"cm")) +
  ylab(expression(hat(beta)[SAR]~-beta))+
  xlab( expression(rho[Y])) + 
  labs(linetype=expression(phi[Y]))


ggsave("FigureA9b.pdf", height=4, width=6, units="in", dpi=600)




